Quantitative comparison of the mean–return-time phase and the stochastic asymptotic phase for noisy oscillators

Seminal work by A. Winfree and J. Guckenheimer showed that a deterministic phase variable can be defined either in terms of Poincaré sections or in terms of the asymptotic (long-time) behaviour of trajectories approaching a stable limit cycle. However, this equivalence between the deterministic notions of phase is broken in the presence of noise. Different notions of phase reduction for a stochastic oscillator can be defined either in terms of mean–return-time sections or as the argument of the slowest decaying complex eigenfunction of the Kolmogorov backwards operator. Although both notions of phase enjoy a solid theoretical foundation, their relationship remains unexplored. Here, we quantitatively compare both notions of stochastic phase. We derive an expression relating both notions of phase and use it to discuss differences (and similarities) between both definitions of stochastic phase for (i) a spiral sink motivated by stochastic models for electroencephalograms, (ii) noisy limit-cycle systems-neuroscience models, and (iii) a stochastic heteroclinic oscillator inspired by a simple motor-control system.

In general, the stochastic case, if it is not a weak perturbation of a deterministic limit-cycle system, requires new phase concepts. Following this insight, the above-mentioned notions of deterministic phase, which are based on Poincaré sections and on the system's asymptotic behaviour, have been generalised to stochastic systems. Whereas Schwabedal and Pikovsky (2013) proposed a notion of phase based on Poincaré sections having a uniform mean-return-time (MRT) property, Thomas and Lindner (2014) introduced an asymptotic phase for stochastic oscillators by means of the argument of the slowest decaying complex eigenfunction of the backward Kolmogorov operator (equivalently, the generator of the Markov process, or the stochastic Koopman operator). However, and in marked contrast to the deterministic case, these two notions of stochastic phase are not equivalent (Cao 2017).
In this paper, we explore the nature of the differences between these two notions of stochastic phase. To this end, we will make use of a recently discovered link between the MRT phase and the Kolmogorov backwards operator (Cao et al. 2020). By exploiting this link, we can calculate both phases using one computational framework; we will use this framework to compare systematically the differences between the MRT phase and the asymptotic phase for a number of stochastic biological systems.
Our paper is organised as follows. In Sect. 2, we review the deterministic phase and isochrons. Next, in Sect. 3 we review both notions of stochastic phase, the MRT and the asymptotic phase. Then, in Sect. 4 we show a procedure relating both expressions which we illustrate by different examples in Sect. 5. We conclude the paper with a discussion of our results. In the "Appendix", we provide details about the numerical methods used in this paper.

The deterministic phase
Consider an autonomous system of ODEṡ whose flow is denoted by φ t (x). Moreover, we assume F(x) is a C 2 vector field having a T -periodic, asymptotically stable, normally hyperbolic limit cycle parameterized by the phase variable θ = 2π t/T Hence, the dynamics of Eq. (1) in Γ can be reduced to a single variable systeṁ As we study attracting limit cycles, any point x ∈ M, where M is the basin of attraction of the limit cycle Γ , will approach Γ as t → ∞ (Hirsch and Pugh 1970). Thus, two points p and q ∈ M will have the same asymptotic phase if This condition extends the notion of phase of oscillation to the basin of attraction M of Γ . Indeed, we can define the function assigning a phase to each point x ∈ M. We can thus define the isochrons as the level sets of ϑ(x), that is which correspond to the leaves of the (strong) stable foliation (that is, the stable manifold M) of Γ (Winfree 1974;Guckenheimer 1975;Winfree 1980). For a normally hyperbolic invariant manifold as M, the phaseless sets correspond to R n \M.

Stochastic phase notions
Next, we review two notions of phase for stochastic systems: the stochastic asymptotic phase and the mean-return-time (MRT) phase. Throughout this Section, we will consider Langevin systems dX dt = f(X) + g(X)ξ(t) where f is an n-dimensional C 2 vector field, g is a C 2 n × k matrix, and ξ is k-dimensional white noise with uncorrelated components ξ i (t)ξ j (t ) = δ(t − t )δ i, j . Moreover, we require the elements g i j (x) in g to be such that the matrix G = 1 2 gg is invertible for all x ∈ R n (see Sect. 4 for further details). For mathematical convenience, we interpret the stochastic differential equation Eq. (7) in the sense of Itô (Gardiner 1985). Schwabedal and Pikovsky (2013) introduced a definition for the phase of a stochastic oscillator in terms of a system of Poincaré sections { MRT (φ), 0 ≤ φ ≤ 2π }, foliating a domain R ⊂ R 2 and possessing a MRT property: a section MRT satisfies the MRT property if for all the points x ∈ MRT the mean return time from x back to MRT , having completed one full rotation, is constant.

The mean-return-time stochastic phase
First defined by Schwabedal and Pikovsky (2013) by means of an algorithmic numerical procedure, the MRT phase was recently related to the solution of a boundary value problem (Cao et al. 2020). As the authors in this paper showed, the MRT sections correspond to the level curves of a function T (x), with appropriate boundary conditions, satisfying the following PDE associated with a first-passage-time problem where L † corresponds to the Kolmogorov backwards operator (the adjoint of the Kolmogorov forward operator L). Both operators read where G = 1 2 gg , and u is an arbitrary C 2 function. Cao et al. (2020) showed that upon imposing a boundary condition amounting to a jump by T (the mean period of the oscillator) across an arbitrary section transverse to the oscillation, the unique solution of Eq. (8), up to an additive constant T 0 , is a version of the so-called MRT function, Hence, the MRT phase Θ(x) satisfies so it evolves in the mean as which is formally analogous to the dynamics for the deterministic phase (see Eq. 3). Thomas and Lindner (2014) defined a notion of stochastic asymptotic phase by means of the eigenfunctions of the Kolmogorov backwards operator. Since the Kolmogorov backwards operator and the stochastic Koopman operator are equivalent (Črnjarić-Žic et al. 2019), the setup in Thomas and Lindner (2014), which we next review, generalises the Koopman approach to obtain the phase of deterministic oscillators to stochastic systems (Mauroy and Mezić 2018;Kato et al. 2021).

The stochastic asymptotic phase
Consider an ensemble of trajectories described by means of the conditional density for s < t. The density evolves following Kolmogorov's equations for L and L † defined in Eq. (9). Assuming the operators L, L † admit a complete biorthogonal eigenfunction expansion with respect to the standard inner product (15) we can write the conditional density as a sum with P 0 suitably normalized, representing the unique stationary probability distribution. The normalization condition in Eq. (16) implies Q 0 ≡ 1. The construction of the stochastic asymptotic phase requires one to assume several properties of the system (7) which Thomas and Lindner termed "robustly oscillatory". First, we require that the nontrivial eigenvalue in Eq. (15) with least negative real part λ 1 = μ + iω is complex, and unique (occurs with algebraic multiplicity one). Thus, we can express its associated right (forward) and left (backward or adjoint) eigenfunctions in polar form as P λ 1 (y) = v(y)e −iφ(y) and Q * λ 1 (x) = u(x)e iψ(x) , where v(y) ≥ 0 and u(x) ≥ 0 are real functions specifying the amplitude of the corresponding eigenfunction.
Second, we require that all other nontrivial eigenvalues λ be significantly more negative, that is, [λ ] < 2μ. This condition guarantees that at sufficiently long times, the sum in Eq. (16) may be written as This asymptotic form means that the density approaches its steady state as a damped focus, with an oscillation period of 2π/ω, and a decaying amplitude with time constant 1/|μ|.
The third assumption is heuristic rather than rigorous: the description as a "stochastic oscillator" will be more appropriate, the larger the quality factor |ω/μ| (Giner-Baldo et al. 2017). That is, provided the oscillation completes sufficiently many rotations before the damping reduces its phase coherence beyond detectability, the system will be "robustly oscillatory". Thus, we require |ω/μ| 1, without specifying an explicit threshold for this quantity.
The dynamics of this focus, capturing the asymptotic oscillatory behaviour of the system, can be obtained in an alternative way: along trajectories X(t), the slowest decaying modes Q * λ 1 (X(t)) evolve in the mean as so they exhibit the same linear focus behaviour as the density ρ(y, t|x, s) when approaching its steady state P 0 . Therefore, we can extract the "stochastic asymptotic phase" ψ(x) from provided u(x) = 0. Analogously to the deterministic case, we will define the points in which a phase cannot be defined as "phaseless sets". The expected value of ψ(x) follows (see "Appendix B" for the complete calculation details) hence, in the limit G → 0, the dynamics for E[ψ(x)] follow the dynamics for the deterministic phase Eq. (3), provided the deterministic system has a well-defined phase (see also "Appendix D" for a discussion about the relationship between ψ(x) and ϑ(x) in the noise vanishing limit). Moreover, if the assumptions under which the uniqueness of solutions of Eq. (12) are met (see next Sect. 4 for a brief review of such conditions), from Eq. (20) it follows that if Following Cao et al. (2020), we assume (without loss of generality) the existence of a parameterisation x = K (α, β) 1

Mathematical relation between the phases
such that the original domain R ⊂ R 2 of system Eq. (7) can be mapped to an annulus via an angular variable α(x) ∈ [0, 2π) and an amplitude-like variable β(x) ∈ [R − , R + ]. Furthermore, we require the noise matrix G in Eq. (9) to be nondegenerate (invertible), so that the L † operator is strongly elliptic (McLean 2000). For the Fokker-Planck equation (9), we impose reflecting boundary conditions at R ± ; for the backward equation (10) we impose adjoint reflecting boundary conditions. For the complete details on the assumptions required for the MRT theory to apply, see Cao et al. (2020).
With these assumptions, we next show how the MRT phase Θ(x) can be represented as the stochastic asymptotic phase ψ(x) plus an additional phase shift Δψ(x). Consider the equality: using Q * λ 1 = u(x)e iψ(x) and the definition of the L † operator, dividing by e iψ(x) , taking the imaginary part and assuming u(x) is a non-vanishing function in R, we find (see "Appendix B" for a complete derivation) Considering the difference between the two phases, which is exactly the definition of MRT phase in Eq. (12). As one can see, since both phase functions, Θ(x) and ψ(x), have a 2π jump, the function Δψ(x) = Θ(x) − ψ(x) will not have a jump. We notice the dependence of L † [Δψ(x)] on the term Δω corresponding to the difference between the frequency ω of the slowest decaying complex eigenmode Q * λ 1 (x) and the MRT frequency 2π/T . Whereas ω can be extracted from the spectra of L † , the MRT period T can be obtained from the stationary probability current J 0 , which is given by where P 0 (x) corresponds to the stationary probability density, satisfying L[P 0 (x)] = 0 and P 0 (x)dx = 1. As shown in Cao et al. (2020), one obtains the mean period by integrating the α-component of the current J 0 along a simple smooth, non-self-intersecting curve, connecting the inner and outer domain boundaries. That is, See "Appendices A.2 and C" for details of the numerical calculation, and an analytical solvable example of Eq. (26), respectively.
In conclusion, solving for and using the first two functions of the eigenvalue problem for L and L † gives us both the stochastic asymptotic phase, the MRT phase, and the difference between the two phases.

Examples
Next, we consider different examples to study how the MRT phase Θ(x) and the stochastic asymptotic phase ψ(x) are related. We will proceed under the assumption that all the models we study satisfy the eigenfunction expansion in Eq. (16), as well as the regularity assumptions given in Cao et al. (2020). For numerical details about computations in this Section, we refer the reader to the numerical "Appendix A".

Spiral sink
We start considering a classical and well-studied stochastic process: a two-dimensional Ornstein-Uhlenbeck process (OUP) in a setting such that the origin becomes a stable sink (Uhlenbeck and Ornstein 1930;Gardiner 1985;Leen et al. 2016;Thomas and Lindner 2019). The general Langevin equation is: where we assume that the two eigenvalues of A are a complex conjugate pair denoted as λ ± = μ ± iω with μ < 0 and ω > 0. We write the matrices A and B as For Eq. (27), the matrix G = 1 2 B B in L † (see Eq. (9)) can be written in the following way Following Thomas and Lindner (2019), we know the asymptotic phase function for Eq. (27) is written as 2 whose expected value follows showing that, as long as there is some noise in the system ( > 0), the term Ω(x) diverges at the origin.
In "Appendix C", we derive the following expression for the mean period T : which, together with Eq. (31) yields that if β D = β c = 0, that is for isotropic noise of the same amplitude, the MRT phase and the stochastic asymptotic phase for the stochastic sink in Eq. (27) are equivalent.
To illustrate the OUP case, we choose coefficients from Powanwe and Longtin (2019), in which a noisy focus is used to model fast gamma-band brain signals (see also Duchet et al. (2020); Spyropoulos et al. (2020) for similar modelling approaches in a neuroscience context). In particular, we take  Figure 1a, b shows both level curves of the stochastic asymptotic phase ψ(x) and the MRT phase Θ(x). As Eq. (31) indicates, the farther we move from the origin, the smaller the term Ω(x) becomes. As a consequence, the difference Δψ(x) is almost negligible far from the origin (panel c) thus causing ψ(x) and Θ(x) to differ just in a small neighbourhood of the origin (panel d). In this case, thanks to the analytical  27). a, b Level curves of the stochastic asymptotic phase ψ(x) and the MRT Θ(x) (y-label shared). c Difference Δψ d Stationary probability distribution (colour coded), with a comparison between some level sets of Ψ (x) (dashed) and Θ(x) (solid). e Comparison of the MRT property (T ≈ 11.07, 2π/ω ≈ 11.13) for ψ(x) (blue) and Θ(x) (orange) for three different MRTisochrons. expression for Δψ (see "Appendix A"), we are able to find the values of the difference near the phaseless set (the origin). As Fig. 1c illustrates, near the origin, the term Δψ(x) alternates between positive and negative values. Panel (e) confirms that the resulting MRT isochrons have the MRT property. While the numerically computed MRT isochrons satisfy the MRT property with high accuracy, the isochrons based on the stochastic asymptotic phase show small but significant deviations from uniformity.

Noisy Wilson-Cowan
Next, we study a noisy version of the Wilson-Cowan (WC) equations, which are widely used to model large-scale neural activity (Wilson and Cowan 1973;Destexhe and Sejnowski 2009;Akam et al. 2012). We adopt the forṁ with S e,i (  Figure 2 displays the results. Panels (a, b) compare the level curves of the stochastic asymptotic phase ψ(x) and the MRT phase Θ(x). In this case, the differences between both level curves are more striking than in the linear focus case. As Fig. 2c shows, and similarly as in the linear focus case, near the phaseless set of the deterministic system, the difference Δψ alternates between positive and negative values. However, the nonlinearities of Eq. (33) cause Δψ to differ from 0 in the bulk of the domain, thus causing differences between the level sets of ψ(x) and Θ(x) (see panel d). Panel (e) demonstrates that the isochrons of the numerically computed MRT phase Θ(x) satisfy the MRT property to within a 1% margin of error, while the isochrons of the stochastic asymptotic phase ψ(x) do not. Moreover, panel (e) also shows that the larger the differences between the level curves of ψ(x) and Θ(x), the larger the deviations of ψ(x) from the MRT property (compare results in panel (e) for I 2π and I 4π Level curves of the stochastic asymptotic phase ψ(x) and the MRT phase Θ(x) (y-label shared). c Phase difference Δψ. d Stationary probability distribution (colour coded), with a comparison between level sets of ψ(x) (dashed) and Θ(x) (solid). e Comparison of the MRT property (T ≈ 6.55, 2π/ω ≈ 6.59) for ψ(x) (blue) and Θ(x) (orange) for three different MRT-isochrons

Noisy Van der Pol oscillator
Next, we study a stochastic version of the Van der Pol equations, which in the absence of noise displays a limit cycle of period T ≈ 6.663. In this case, we set the noise in each component 1. Due to their similarity with the FitzHugh-Nagumo model (FitzHugh 1961;Nagumo et al. 1962), these equations are widely used in neuroscience as a useful reduction in the Hodgkin-Huxley neuron model (Izhikevich 2007). At a macroscopic level, they are also used to describe successfully the dynamics of epileptic tissue (Proix et al. 2017; Pérez-Cervera and Hlinka 2021).
We illustrate results for this oscillator in Fig. 3. As panels (a, b) illustrate, differences between the stochastic asymptotic phase ψ(x) and the MRT phase Θ(x) are very small and they are restricted to a neighbourhood of the origin. Indeed, the structure of the discrepancies is similar to the differences for the OUP. That is, they alternate between negative and positive values (compare panel (c) in Figs. 1, 3). We observe (panel d) that both phases are almost identical near the maxima of the stationary distribution P 0 . Together with the near equivalence of both periods 2π/ω ≈ 6.59, T ≈ 6.55, their similarity causes both phases to be nearly indistinguishable in this case. Indeed, as panel (e) shows, the MRT property is satisfied very accurately for both phases, except for ψ(x) in points very near the origin (where phase discrepancies Δψ(x) are larger).

Noisy heteroclinic oscillator
We complete our comparison between the stochastic asymptotic phase ψ(x) and the MRT phase Θ(x) by studying a noisy heteroclinic oscillator. The specific form of the deterministic heteroclinic system we consider was introduced in Hirsch et al. (2012), chapter 10, and was adapted to a biological context as a conceptual model for a central pattern generator control mechanism based on a dynamical architecture alternative to the standard limit cycle architecture (Shaw et al. 2012(Shaw et al. , 2015Lyttle et al. 2017;Park et al. 2018). We study the form given in (Giner-Baldo et al. 2017), namelẏ with α = 0.1, D = 0.01125 and reflecting boundary conditions on the domain −π/2 ≤ {X, Y} ≤ π/2. Without noise, Eq. (35) has an attracting heteroclinic cycle, consisting of a closed loop of trajectories connecting a sequence of saddle equilibria which is capable of sustaining robust oscillations in the presence of noise (Thomas and Lindner 2014). We study this oscillator for lower (D = 0.01125) and higher (D = 0.1) level of noise and present results for each case in Figs. 4 and 5, respectively. Panels (a, b) compare the stochastic asymptotic phase ψ(x) and the MRT phase Θ(x). In contrast to the previous cases, the structure of the phase difference Δψ(x) appears to be rotationally symmetric, to a good approximation, in a neighbourhood of the centre (see both c panels). These differences in phase lead to the principal discrepancies between the level curves of both phases appearing near the origin (panel d). Indeed, as the MRT property check in panel (e) shows, whereas the MRT property holds quite well for points computed using the MRT phase Θ(x), it does not hold in general for the stochastic asymptotic phase ψ(x). However, for lower noise, the MRT property holds for points away from the origin in which the level curves for both phases coincide and the stationary probability is concentrated. As noise increases, the mean return time of the stochastic asymptotic phase becomes increasingly positiondependent close to the domain border.

Discussion
In this paper, we have derived a framework to compute simultaneously the MRT phase and the asymptotic phase of a stochastic oscillator. Our results build on prior work by Schwabedal and Pikovsky (2013) and Thomas and Lindner (2014) defining the notions of MRT phase Θ(x) and stochastic asymptotic phase ψ(x), respectively. While initially defined on an algorithmic basis, the MRT phase was recast by Cao et al. (2020) as the solution of a PDE with jump-periodic boundary conditions. As a result of Cao et al. (2020), a relationship between the Kolmogorov backwards L † operator and the MRT phase was derived (see Eq. (12)). Since the stochastic asymptotic phase was already defined as the argument of the slowest decaying eigenfunction of L † , in this work we developed the link between both phases and the Kolmogorov backwards operator to obtain an expression for the difference between the two phases. That is, The computation of the difference Δψ(x) allowed us to compare both phases in different dynamical scenarios: we have considered two examples of noise induced oscillations (a spiral sink and an heteroclinic oscillator) and  20) is not zero. As Cao (2017) observed, for a planar oscillator with isotropic noise, the condition Ω(x) = 0 is satisfied if the eigenfunction Q * λ 1 is a complex analytic function of its arguments (in the sense of complex variables theory). But the practical significance of the difference Δψ(x) has not been systematically explored before now.  articulate the distinction between pathwise and ensemble descriptions of a stochastic process. For example, a general diffusion process may be described either in terms of the Itô stochastic differential equatioṅ

Two perspectives on stochastic oscillators
which describes the evolution of a single trajectory along a sample path, or in terms of the Fokker-Planck equation which describes the evolution of the density ρ(x, t) = 1 |dx| Pr [X(t) ∈ [x, x + dx)] of an ensemble of trajectories (Øksendal 2003). Similarly, a discrete chemical reaction process comprising M reactions with stoichiometry vectors ζ k and hazard functions λ k , driven by independent Poisson processes {Y k (t)} M k=1 , has a pathwise description of the form as well as an evolution equation (the so-called chemical master equation (Higham 2008)) of the form where p(x, t) is the probability that the state X(t) is exactly x at time t (Anderson and Kurtz 2015). The transit time for a stochastic oscillator to reach a Poincaré section from a starting point on that section, having completed one full rotation, is a stopping time (Karatzas and Shreve 2012), thus, a random variable that arises from the individual sample path. The "mean-return-time" function T (x) (Cao et al. 2020) is defined from the ensemble average of this quantity. Importantly, the MRT property describes the behaviour of trajectories over a finite time horizon, namely looking roughly one period into the future.
In contrast, one defines the stochastic asymptotic phase ψ(x) (Thomas and Lindner 2014) in terms of the longtime statistical behaviour of an ensemble of trajectories, as captured by the biorthogonal eigenfunction expansion Eq. (15) of the forward and backward operators. Thus, while the equations satisfied by both the MRT function, namely L † [T ] = −1, and the stochastic asymptotic phase eigenfunction, namely L † [Q] = λQ, involve the backward Kolmogorov operator, we see the MRT as related to a "pathwise" description over a finite time horizon, and the asymptotic phase as related to the "ensemble" description of the process at long times.
We speculate that this distinction may turn out to play a role in choosing which notion of phase applies more naturally to specific problems, such as synchronization of coupled stochastic oscillators (long-time behaviour), or "phase response" of a stochastic oscillator to a single kick (short-time behaviour). Fortunately, as we have seen above, for many biological examples the quantitative difference between the two types of phase is small. Their quantitative similarity thus provides investigators a degree of flexibility in working with whichever notion of phase is conceptually best suited to a given problem.

Noise amplitude
For systems with an underlying limit cycle, we observe empirically that both phases appear to coincide with the deterministic phase as the level of noise approaches zero. For the MRT phase, its convergence to the deterministic phase in the case of vanishing noise has been investigated by Cao et al. (2020) §2.4, who established convergence under additional regularity assumptions. For the stochastic asymptotic phase, its convergence to the deterministic phase for vanishing noise was anticipated by Thomas and Lindner (2014), see also the discussion in terms of the Koopman operator by Kato et al. (2021) and derived in "Appendix D" in this manuscript using the same additional regularity assumptions as in Cao et al. (2020). In line with these observations, although the intrinsic differences between both phases depend on each particular system, we have found the differences between the MRT phase and the stochastic asymptotic phase to grow as the noise is increased. It is nevertheless interesting to observe that, in every case we have explored, the differences between both phase level curves were restricted to areas where the stationary density was low. Therefore, once differences in mean period were accounted for, both phases were practically indistinguishable when describing single trajectories, differing only in how well they satisfy the MRT property, which is a property of the ensemble.
For a deterministic limit cycle (LC) system, in which both phase perspectives coincide, the function e iψ(x) is an eigenfunction of L † with eigenvalue λ 1 = iω. As soon as some noise is introduced in the system (G = 0), the eigenvalue λ 1 associated with the slowest decaying eigenfunction becomes complex instead of purely imaginary, that is λ 1 = μ + iω (with μ < 0). As a consequence, for G = 0, the information about the initial phase is dissipated as time progresses. By contrast, we can think of the MRT function as a function containing information about the oscillatory system which does not vanish as t → ∞. Therefore, for systems having an underlying LC, one can expect that the more robustly oscillatory the system (i.e. |μ| ω), the more similar ψ(x) and Θ(x) become. This interpretation agrees for the LC systems we studied: the Wilson-Cowan (WC) equations and Van der Pol (VdP) model. For the WC equations, we observed larger discrepancies between ψ(x) and Θ(x) than in the VdP model. These differences cause a loss of the MRT property for the stochastic phase, which is consistent with the magnitude of |μ/ω| for both cases: |μ/ω| ≈ 0.415 for the WC equations and |μ/ω| ≈ 0.054 the VdP model, respectively (see Table 1 and Fig. 6 at "Appendix A").
Despite not having a limit cycle, we observe that for the noisy heteroclinic oscillator (HO), the system is less robustly oscillatory as the noise increases. More precisely, |μ/ω| ≈ 0.115 for D = 0.01125, and |μ/ω| ≈ 0.269 for D = 0.1, respectively (see Table 1 and Fig. 6). Indeed, our computations for the HO suggest that both μ and ω tend to 0 as G → 0, which can be interpreted as the system approaching an infinite period closed loop as G → 0. Hence, the results for the HO can be considered a particular case of the LC case, thus supporting the interpretation of the role of |μ/ω| in the loss of the MRT property for the stochastic asymptotic phase.
The focus case requires a different interpretation. Unlike the LC or HO cases, stable focus systems have no closed loop structure in the absence of noise. Indeed, in the focus case μ, approaches the negative nonzero real part of the stable focus eigenvalue as G → 0. As a consequence, the addition of noise may, in general, increases or decreases |μ/ω|. In this paper, we considered a widely used linear model: the Ornstein-Uhlenbeck process (OUP). For this model, we provided a new formula T and an initial seed for computing Δψ(x) accurately even near the origin, thus facilitating the computation of its MRT function. However, the OUP system we studied has the particular property of not changing |μ/ω| as the noise increases. Studying how the change of μ/ω affects the phase dynamics and other important features of nonlinear stochastic foci, such as the amplitude of the stochastic oscillator (Pérez-Cervera et al. 2021), appears as an interesting topic for further research.

Future perspectives
In developing our numerical examples, we have proceeded under the assumption that the robustly oscillatory criteria are met and that each system studied has a complete biorthogonal eigenfunction expansion. Whether this is rigorously true or not in specific cases is a question of functional analysis that goes beyond the scope of this paper. However, although our theoretical development assumes the expansion, it appears that practically, all that is really required is that the low-lying spectral elements (eigenmodes with eigenvalues having relatively small real and imaginary parts) exist and are discrete.
From a numerical perspective, the methodology introduced in this paper extends the numerical procedure introduced in Cao et al. (2020) which assumes the location of the phaseless set to be known a priori (see Section A.4 in the "Appendix" for further discussion of this point). However, whereas our procedure is restricted to systems in which the set of SDEs describing the system is known, and the noise is temporally uncorrelated and Gaussian, the procedure presented in Schwabedal and Pikovsky (2013) (based on Monte Carlo simulations) applies to a wide range of systems, noise types and also to data. Nevertheless, due to the equivalence between the Kolmogorov backwards operator and the stochastic Koopman operator (Črnjarić-Žic et al. 2019), we expect that the computation of the eigenfunctions of interest from data via dynamical decomposition methods (Schmid 2010;Budišić et al. 2012;Proctor et al. 2016;Brunton and Kutz 2019;Mauroy et al. 2020), combined with the theoretical expression we gave in Eq. (24), may provide an alternative way of obtaining the MRT phase from data.

A Numerical details
In this Appendix, we detail the numerical methodology leading to the results in this paper. Results in this Appendix follow from previous work in Cao et al. (2020);.

A.1 Diagonalization of the L † operator
We construct the L † operator as follows.
Given a Langevin equation as in Eq. (7), we restrict its phase space to a rectangular domain whose size is chosen large enough so that the probability for individual trajectories X(t) to reach the boundaries is very low. For the special case of the heteroclinic oscillator, boundaries are given by the nature of the system. Then, we discretize the domain D in N and M points such that Δx = (x + − x − )/N and Δy = (y + − y − )/M. Then, we build L † by using a standard centred finite difference scheme, except at the borders of the domain. For the heteroclinic oscillator, we implemented adjoint reflecting boundary conditions at the borders of the domain Gardiner (1985). In contrast, for the unbounded systems, since there is no natural border, we substitute the centred finite difference scheme by a forward (or backward) finite difference scheme over a bounded domain. Using adjoint reflecting boundary conditions for these systems yielded numerically very similar results. By diagonalizing the resulting matrix, we obtain the eigenvalues and the associated eigenfunctions of L † . We remark we are not interested in the complete spectrum of L † but on the small part of it (see Fig. 6). As we review at Sect. 3, we just need to consider the eigenvalue associated with the slowest decaying complex eigenfunction Q * λ 1 (x) to obtain the functions u(x) and ψ(x).

A.2 Computing the MRT period T
To compute the MRT period T , we build the Kolmogorov forward operator L in Eq. (9) following the same numerical procedure as in A.1 underlying the construction of L † . We obtain P 0 as the eigenfunction of L with null eigenvalue. Then, as we explained in Sect. 4, we compute the integral . The stochastic asymptotic phase ψ(x) can be obtained from the argument of the eigenfunctions Q * λ1 (x) associated with the eigenvalues the smallest non-negative purely complex eigenvalues λ ± (λ + = λ 1 ). In addition, the eigenfunction Σ(x) associated with the smallest non-negative purely real eigenvalue λ Floq corresponds to the stochastic isostables. Right panel: level curves of the isostable function Σ(x), ten trajectories and the stochastic limit cycle Σ 0 ≡ {x|Σ(x) = 0} (black curve) in Eq. (26) to obtain T . More precisely, if we denote the phaseless point as (x,ȳ), we determine T by integrating the y component of the stationary probability current J 0,y along the line joiningx and Of course, other choices may be possible (e.g. integrating the x component along the y axis).

A.3 Computation of the phase offset 1Ã
In order to obtain the phase difference Δψ, we interpret the system as a linear system of the form Ax = B. Since we have already build the L † operator, and we already know ω and T , we just need to build the terms Ω(x) and Δω. The computation of the derivatives of ln(u(x)) and ψ(x) in Ω(x) is done using a finite difference scheme. However, solving above Eq. (38) requires to deal with two sources of numerical instability. On the one hand, the derivatives of ln(u(x)) and ψ(x) diverge at the phaseless set. For this reason, the computed numerical values of Ω(x)+Δω near the phaseless set are not accurate. On the other hand, the L † operator is represented by a sparse matrix, thus making the computation of its inverse numerically unstable. To overcome these numerical issues, we remove the values of Ω(x) + Δω falling inside a small radius r min see (Table 1) around the phaseless set and solve Eq. (38) by means of a least squares iterative method to approximate the solution. This procedure finds a solution Δψ(x) of Eq. (38) having a very small error (for the studied cases in this paper max error values were O(10 −3 )). Then, we obtain an estimate for Δψ(x) for the points of the grid inside r min by using extrapolation routines.
For the noisy spiral sink in Eq. (27), we found with γ = −β c μ + β D ω and α = β c ω + β D μ to be a very good initial seed since it solves Eq. (38) with an O( 2 ) error. Thanks to this initial seed, it was not necessary to remove any point from the least square iterative solver (that is to say r min = 0).

Fig. 7
Comparison for different topologies between the MRT sections for the noisy heteroclinic oscillator with isotropic noise (D = 0.01125). The blue-yellow colour scheme corresponds to the results published in Cao et al. (2020) (with a small hole around the origin, corresponding to an annulus). We superimposed the results of the present paper (see Fig. 4c) without a hole (hence, a topological grid). Both topologies yield numerically indistinguishable isochrons a sufficiently high degree of symmetry, the location of the phaseless set may be clear a priori. For instance, in the heteroclinic oscillator, which has the symmetry of the square, the phaseless point should be at the centre of the square. However, in other systems, such as the Wilson-Cowan system, the location may not be known a priori. This lacuna prevents one from constructing an annular domain numerically that is guaranteed to exclude the phaseless point. Fortunately, we have observed that in systems for which the phaseless set's location is known, there is no practical difference in the structure of the isochrons produced by solving the MRT equation with a small central exclusion with reflecting boundary conditions, and a construction with a grid that covers the entire domain without implementing the central excluded region (see Fig. 7). This apparent robustness of the numerical procedure means that in a nonsymmetric system such as the Wilson-Cowan example, we may construct our numerical implementation without the inner annular exclusion. The resulting isochrons will converge and conflict in a small region on the scale of a single grid spacing. This singularity localizes the phaseless point to within the accuracy of a grid spacing.
For these reasons, despite the topological differences, we prefer to use a complete rectangular grid instead of an annulus for numerical implementation. This method leads to a simpler implementation which nevertheless still yields numerically indistinguishable results, with respect to the ones which would be obtained by explicitly adding a hole around the phaseless set. On a theoretical level, there are four separate kinds of isochrons that one may consider: the eigenvalue λ and eigenfunction Q = ue iψ for the simply connected domain, and for an annular domain; and the MRT phase θ for the annular domain and the MRT phase obtained numerically for a simply connected domain. In this paper, we investigated analytically the relationship between ψ and θ on a general domain and investigated numerically the two phases on the simply connected domain. The MRT isochrons obtained using the (theoretically questionable) simply connected domain nevertheless satisfied the MRT property: the mean time to return to the isochron was independent of the starting position along the isochron.
Thus, the topological discrepancy between the annular and simply connected domain does not seem to be significant numerically, at least for the examples we considered in this paper. It is an interesting open question to ask how λ and Q compare on the annular versus the simply connected domain. What impact does the size and location of the annular exclusion have on λ and Q? Moreover, how do λ and Q change if the hole is located at a point not overlapping the phaseless set of the simply connected domain? In cases for which we do not know a priori the location of the phaseless set, these questions may become highly relevant.

A.5 Choosing the zero-th phase through the stochastic isostables
Like deterministic phase variables, the MRT phase Θ and the stochastic asymptotic phase ψ are defined up to an arbitrary additive constant. In deterministic limit cycle systems, the phase is often chosen to be zero at the maximum of some variable of interest, such as the voltage of a spiking neuron. In order to compare Θ and ψ, we use a recently introduced generalization of the isostable coordinate adapted to stochastic systems . Briefly, just as the slow-est decaying complex eigenmode Q * λ 1 (x) allows one to define a stochastic phase ψ(x), the slowest decaying purely real eigenmode allows one to define a stochastic amplitude since it accounts for the slowest mode describing pure contraction without an associated oscillation. Hence, the level curves of such a slowest decaying purely real eigenmode-which we denote as Σ(x)-correspond to the stochastic isostables (see Fig. 6). For this reason, and following the usual approach from deterministic systems, in this paper we used the maximum of the "stochastic limit cycle", corresponding to the 0-level isostable Σ 0 ≡ {x | Σ(x) = 0} in the x direction, to set the zero phase point for both phases Θ(x) and ψ(x). In this way, we are able to provide a consistent basis for comparison throughout the paper.

A.6 Computation of the MRT property
Once Δψ is computed, we obtain the MRT phase as Θ(x) = ψ(x) + Δψ(x). To check if the computed function Θ(x) satisfies the MRT property, we do the following for each point x 0 : first, we interpolate Θ(x) in the whole domain D in Eq. (36) by using a 2D spline method. Then, we integrate the system of interest for a time large enough to assure that the first return occurs with overwhelming probability. In practice, we integrate for a time 3T by means of a Euler-Heun scheme with a time step O(10 −3 ). By using the interpolated grid, we can obtain a description of the trajectory X(t) in terms of the phase θ(t) = Θ(X(t)). Hence, we can check at which time we first cross Θ(x * ) = Θ(x 0 )+2π . We repeat the procedure averaging the return time over 10 5 realisations. To check the MRT property for the stochastic asymptotic phase ψ(x), we repeat the previous procedure just substituting Θ(x) by ψ(x). The computed standard error of the MRT results in panel (e) of Figs. 1, 2, 3, 4 and 5 was found to be less than 0.1.

B Derivation of d dt E[Ã(x)] (Eq. 20) and L † [Ã(t)] (Eq. 23)
In this section, we give details of the derivation of Eqs. (20) and (23). We thank the anonymous reviewer who provided the following elegant derivation.
For completeness, we also state whose significance remains to be explored elsewhere.
Then, the probability current J can be written as so the mean period T can be computed as which yields from which we obtain T in Eq. (32).

D Relation to the asymptotic phase for deterministic systems
In this section, we explore the relationship between the mean return time (MRT) phase Θ(x), the stochastic asymptotic phase ψ(x), and the deterministic phase ϑ(x). Under certain regularity assumptions, Cao et al. (2020) §2.4 established convergence of Θ(x) to ϑ(x) in the limit of vanishing noise.
Here, we present a similar argument as in Cao et al. (2020) §2.4, to establish that if suitable regularity conditions are satisfied, then in the limit of small noise the stochastic asymptotic phase ψ(x) likewise converges to ϑ(x).
We start by taking the time derivative of the deterministic phase function ϑ(x), defined in Eq. (5), which holds ∀x in the basin of attraction of the limit cycle. Consider a family of stochastic differential equations as in Eq. (7), but with the noise scaled by a small parameter √ . That is, where ξ is a vector with components comprising independent delta-correlated white noise. Moreover, we also consider the corresponding family of solutions of Eq. (20), that is We now make the regularity assumption that as → 0, ψ (x) converges uniformly on compact subsets of the domain to a C 2 function ψ 0 (x). 3 As for any , ψ (x) is defined up to an additive constant, we consider convergence in the sense that for arbitrary nonzero vectors v ∈ R 2 , v (∇ψ 0 (x) − ∇ψ (x)) → 0 as → 0, for all x in the domain. Fixing x and setting v = f(x), we see for each x as → 0; here, we have used Eq. (53). Consequently, if ψ (x) converges to a well-behaved function ψ 0 (x) in this way, it must satisfy Comparing Eqs. (51) and (56), evidently if the deterministic system has a stable limit cycle, then the function ψ 0 (x) must correspond with the deterministic phase ϑ(x) through the linear relation for an arbitrary constant ϑ 0 . Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.